% Copyright (C) 2009,2010,2011,2012  Marco Restelli
%
% This file is part of:
%   LDGH -- Local Hybridizable Discontinuous Galerkin toolkit
%
% LDGH is free software: you can redistribute it and/or modify it
% under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% LDGH is distributed in the hope that it will be useful, but WITHOUT
% ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
% or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
% License for more details.
%
% You should have received a copy of the GNU General Public License
% along with LDGH. If not, see <http://www.gnu.org/licenses/>.
%
% author: Marco Restelli                   <marco.restelli@gmail.com>


function [Xi,Ti,Ui,cell_type] = interpolate_grad(grid,base,k,U,varname)
% [Xi,Ti,Ui,cell_type] = interpolate_grad(grid,base,k,U,varname)
%
% This function is analogous to interpolate, except for the fact
% that the gradient is computed.
%  varname: can be {'phi'}
% The results are as follows:
%  Xi: nodal coordinates
%  Ti: connectivity of the reconstructed field
%  Ui: interpolated values
%  cell_type: paraview code for the elements used in the reconstructed
%             field
%
% See interpolate for further details.

 % 1) Get the nodes
 locmeshes = load('localLagmeshes.octxt');
 if(k>size(locmeshes.meshes,2))
   error(['Degree ',num2str(k),' not supported.']);
 endif
 switch varname

  case {"phi"}
   locmesh = locmeshes.meshes(base.me.d,k);

  otherwise
   error (["Unknown varname value ",varname]);
 endswitch

 p = locmesh.p;
 t = locmesh.t;
 cell_type = locmesh.cell_type;

 % 2) Evaluate the basis on the local mesh
 np_tot = size(p,2); % number of points per old cell
 np     = size(t,1); % number of points per new cell
 nt     = size(t,2); % number of new cells per old cell
 switch varname

  case "phi"
   gPHI = zeros(grid.d,base.pk,np_tot);
   for i=1:base.pk
     for id=1:grid.d
       gPHI(id,i,:) = ev_pol(base.gradp_s{id,i},p);
     end
   end

 endswitch

 % 3) Build nodes, connectivity and nodal values
 idx = ones(1,np);
 
 switch varname

  case {"phi"}
   Xi = zeros(grid.d,np,nt*grid.ne);
   Ti = zeros(       np,nt*grid.ne);
   switch varname
    case "phi"
     Ui = zeros(grid.d,np,nt*grid.ne);
   endswitch
   for ie=1:grid.ne
     for it=1:nt

       i = (ie-1)*nt + it;
       Xi(:,:,i) = grid.e(ie).b * p(:,t(:,it)) + grid.e(ie).x0(:,idx);

       Ti(:,i) = (ie-1)*(nt*np) + (it-1)*np + [1:np];

       switch varname
        case "phi"
	 for id=1:grid.d
           Ui(id,:,i) = (permute(gPHI(id,:,t(:,it)),[2,3,1])' * U(:,ie))';
	 end
	 % Now the coordinate transformation
	 Ui(:,:,i) = grid.e(ie).bi' * Ui(:,:,i);
       endswitch
     end
   end

 endswitch

return


